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PACS 05 . 45 . Jn - High dimensional chaos 

PACS 05 . 45 . Xt - Synchronization 

PACS 31.30-Jh- Long-range interactions 

PACS 05 . 70 . Fh - Phase transitions in statistical mechanics and thermodynamics 

Abstract. - We study the chaotic behavior of the synchronization phase transition in the Ku- 
ramoto model. We discuss the relationship with analogous features found in the Hamiltonian 
Mean Field (HMF) model. Our numerical results support the connection between the two mod- 
els, which can be considered as limiting cases (dissipative and conservative, respectively) of a 
more general dynamical system of damped-driven coupled pendula. We also show that, in the 
Kuramoto model, the shape of the phase transition and the largest Lyapunov exponent behavior 
are strongly dependent on the distribution of the natural frequencies. 



X 



, Introduction. — Long-range interacting systems 
have been intensively studied in the last years and new 
methodologies have been developed in the attempt to un- 
derstanding their intriguing features. One of the most 
promising directions is the combination of statistical me- 
chanics tools and methods adopted in dynamical systems 
[X\. In particular, phase transitions have been extensively 
explored in both conservative and dissipative long-range 
^stems. The Hamiltonian Mean Field (HMF) model [2] 
and the Kuramoto model [3-5] represent two paradigmatic 
tioy models, the former conservative and the latter dissi- 
pative, for many real systems with long-range forces and 
Ijave several applications. Both models share the same 
order parameter and display a spontaneous phase transi- 
tion from an homogeneous/incoherent phase to a magne- 
tized/synchronized one. 

In [6] we already observed that HMF and Kuramoto mod- 
els can be considered as limiting cases (respectively con- 
servative and overdamped) of a more general model of 
driven-damped coupled incrtial oscillators. In this paper 
we present new numerical results which support a common 
scenario for the two models. More precisely, first we dis- 
cuss the well known equilibrium features of the second or- 
der phase transition in the HMF model, then we study the 
stationary asymptotic behavior of the Kuramoto model as 
a function of the coupling strength. On one hand, through 
new numerical simulations of large size systems, we con- 



firm that, as also pointed out by other authors [7] , the 
shape of the dynamical phase transition in the Kuramoto 
model changes from a continuous to an abrupt one, de- 
pending on the distribution of the natural frequencies of 
the oscillators. On the other hand, and this is our main 
result, we clearly show that, as for the HMF model, the 
largest Lyapunov exponent (LLE) of the Kuramoto model 
exhibits a peak just around the critical value, confirming 
the generality of this microscopic signature for a phase 
transition. Chaotic behavior in the Kuramoto model was 
discussed previously in ref. [8]. However those authors 
compute the entire spectrum of the Lyapunov exponents 
only for small sizes and only for one kind of natural fre- 
quencies distribution, without discussing the strong de- 
pendence of the chaotic behavior on the shape of that dis- 
tribution and its persistence in the thcrmodynamical limit. 
Finally, by tuning the width of the natural frequencies dis- 
tribution, we show how the phase transition changes conti- 
nously from 2nd-order-like behavior towards a Ist-order- 
like one and we draw a complete synchronization phase 
diagram. As far as we know, these results are reported for 
the first time and we think that they could provide with 
new insights for the study of dynamical phase transitions 
in systems displaying collective synchronization. 

Phase transition and chaos in the HMF model. 

The Hamiltonian Mean Field model describes the dynam- 
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ics of A^ classical spins or rotators, characterized by the an- 
gles ^i £ [—TT, 7r[ and the coniugate momenta Pi €] — oo, cxd[, 
which can also be represented as particles moving on the 
unit circle. In its ferromagnetic version the Hamiltonian 
of the model is given by: 



H = K + V 



yz^ + ^y [1 



,)], (1) 



i,i = l 



where i — 1, ..., A^ and the mass m is usually set to 1. 
The potential term of Eq.l reveals the mean field nature 
of the model, since each rotator can interact with all the 
others. Such a nature becomes more evident if we de- 
fine as order parameter the magnetization M — A/e"^ = 



^■Ef 



JOi 



where M and 6 are the modulus and the 



N ^3 = 1 ' 

global phase. Within this assumption the Hamilton equa- 
tions of motion can be written 



1 " 
N ^ 



,) = A/sin( 



,), i-1, ...,A, (2) 



which correspond to the equations of single pcndula in a 
mean field potential. We note also [6] that Eq.2 can be 
regarded as the conservative limit of the following mean 
field equation describing a system of driven and damped 
pendula (with unit mass): 



6*, + 561; + CA/ sin (61, 



1, ...,A, (3) 



provided that the coupling C = 1, the damping coefficient 
5 = and the torque term F = 0. 

The equilibrium solution of the HMF model can be derived 
in both the canonical and microcanonical ensembles [2]. 
It gives the exact expression of the so-called caloric curve, 
i.e. the dependence of the energy density U ~ H/N on the 
temperature T: U = i^ + \{l- M^), being f3 = 1/T, and 
predicts a second-order phase transition from a low-energy 
condensed phase (with M > 0) to an high-energy homo- 
geneous one (with M = 0) at the critical temperature 
Tc ~ 1/2 (corresponding to the critical energy density 
C/ = 4/3). 

The microcanonical simulations at equilibrium confirm 
these predictions and also allow to get some information 
about the microscopic dynamics of the system [2,9] . It is 
well known that in classical many-particle systems macro- 
scopic collective behavior can coexist with chaos at the 
microscopic level. This feature is particularly evident near 
a phase transition, where chaotic dynamics can induce 
non trivial time dependence in macroscopic quantities. In 
these cases it is worthwhile to study the Largest Lyapunov 
exponent (LLE), which gives a sufficient condition for 
chaotic instability by measuring the asymptotic rate of ex- 
ponential growth of vectors in tangent space. For this pur- 
pose one has to consider the limit: A = limt_ 

where d(t) = ^/Y~m¥^^W?) 

tance calculated from the infitesimal displacements at time 



t ^^^ d(0) ' 

is the Euclidean dis- 
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Fig. 1: (Lower panel) The modulus of the order parameter M of 
the HMF model is plotted as a function of the energy density U 
for various system sizes at equilibrium: A^ = 1000, 5000, 10000; 
(Upper panel) Numerical calculation of the equilibrium largest 
Lyapunov exponent (LLE) as a function of U for the same 
system sizes. In both the panels each dot represents an average 
over 10 runs. See text. 



t. Then, in order to obtain the time evolution of d(t\ 
one must integrate along the reference orbit the linearized 
equations of motion following for example the procedure 
of ref. [10]. 

In the upper panel of Fig.l we plot the LLE as a function 
of the energy density for increasing system sizes A^, while 
in the lower panel the correspondent magnetization curve, 
exibiting the typical shape of a continuous phase transi- 
tion, is reported for comparison. An average over 10 real- 
izations at equilibrium has been considered for each point. 
As expected, in both the limits of small and large ener- 
gies, where the system becomes integrable, the LLE goes 
to zero. On the other hand, just before the critical energy, 
LLE exhibits a peak which persists and becomes broader 
increasing the size A" (see also [9]). In particular, it has 
been already shown (see Fig. 16 in [2]) that the LLE is pos- 
itive and N-independent just below the transition, while 
it goes to zero above it (as A'"^'"^) and also for very small 
energy densities. Such a behavior at equilibrium is strik- 
ingly correlated to the kinetic energy fluctuations [2,9] and 
it is also in agreement with a theoretical formula relating 
the LLE with other dynamical quantities [11], see [12,13] 
for the general theory. In the next section we will show 
that analogous features can be found also in an apparently 
different context, as that one of the Kuramoto model. 

Phase transition and chaos in the Kuramoto 
modeL — The Kuramoto model [3-5] is considered one 
of the simplest models exhibiting spontaneous collective 
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Fig. 2: (Lower panel) The asymptotic order parameter r of the 
Kuramoto model is plotted as a function of the ratio K/Kc 
for several system sizes [N = 20000, 30000, 40000 and 50000) 
and for a Gaussian distribution g{uj) of the natural frequencies. 
Decreasing values of KjKc are used in order to compare the 
data with the HMF ones. In this case Kc = 1.59617... A 2nd- 
order-like dynamical transition from the homogeneous phase to 
a synchronized one, similar to that one observed in the HMF 
model (Fig.l), is clearly visible. (Upper panel) The Largest 
Lyapunov Exponent (LLE) is plotted as function of K/Kc- 
A well defined peak around the phase transition indicates a 
microscopic chaotic activity which is maximal at the critical 
point 
. In both panels each dot represents an average over 10 
realizations. See text. 



Fig. 3: (Lower panel) The asymptotic order parameter r of the 
Kuramoto model is plotted as a function of the ratio K/Kc for 
several system sizes (iV = 20000, 30000, 40000 and 50000) and 
for a uniform distribution g{ijj) of the natural frequencies. The 
critical coupling is in this case Kc = S/tt. A Ist-order-like 
dynamical transition from the homogeneous phase to a syn- 
chronized one is clearly visible. This is more evident in the 
inset where a "symmetric" uniform distribution was used, see 
text. (Upper panel) The Largest Lyapunov Exponent (LLE) 
is plotted as function of K/Kc- A sharp peak around the 
phase transition indicates microscopic chaotic activity. Also 
in this case we show in the inset the results for the "symmet- 
ric" frequency distribution, see text, for which the transition is 
sharper. In all panels each dot represents an average over 10 
runs. See text. 



synchronization. It describes a large population of cou- 
pled limit-cycle oscillators, each one characterized by a 
phase Oi and a natural frequency uji, whose dynamics is 
given by: 

^' = ^' + ]7Z^«' 



accordance vifith the expression Kc 



3(oy 



sm 



% - 00, 



(4) 



where K > is the coupling strcnght and i = l,...N. 
The natural frequencies are time- independent and arc ran- 
domly chosen from a symmetric, unimodal distribution 
g{uj). We will consider here only uniform and Gaussian 
g{uj) distributions. As in the case of HMF model, one 
can immagine the oscillators as particles moving on the 
unit circle. For small values of K, the oscillators act as 
if they were uncoupled and each oscillator tends to run 
independently and incoherently with its own frequency. 
Instead, when K exceeds a certain threshold Kc, the cou- 
pling tends to synchronize each oscillator with all the oth- 
ers and the system exhibits a spontaneous transition from 
the previous incoherent state to a synchronized one, where 
all the oscillators rotate at the same frequency fl (a value 
which corresponds to the average frequency of the sys- 
tem, preserved by the dynamics). As shown by Kuramoto 
itself [3], the critical value of the coupling depends only 
on the central value g{Lu = 0) of the distribution g{Lo) in 



The order parameter of the Kuramoto model is perfectly 
equivalent to the magnetizaton in the HMF model and it 
is given by r = re'^^ = ^ ^7=1 ^'^^ , where (j) is, again, 
the average global phase corresponding to the centroid of 
the phases of the oscillators and the modulus < r < 1 
represents the degree of synchronization of the population. 
In terms of the variables r and (j), eq. (4) can be rewritten 
as: 

0,=:uj^ + Krsin{(l)-0,), i = l, ...,7V. (5) 

where, as happened also for the HMF model, the mean 
field character of the system becomes obvious. For a given 
value of K, as the population becomes more coherent, r 
grows and the effective coupling Kr increases. In this 
regime of partial synchronization, as predicted by the so- 
lutions of eq. (5) , two kinds of oscillators coexist depend- 
ing on the size of \uJi\ relative to Kr: (i) oscillators with 
\LUi\ > Kr, called drifting — oscillators, that run incoher- 
ently around the unit circle; (ii) oscillators with \iOi\ < Kr, 
called locked — oscillators, that are trapped in a rotating 
cluster. The dynamic interplay between these two kinds 
of oscillators is probably at the root of the microscopic 
chaotic behavior which, as we will show, characterizes the 
regime of partial synchronization. On the other hand, 
when the effective coupling Kr becomes strong enough, all 
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the oscillators rotate in the same cluster at the frequency 
n and any fingerprints of chaos disappears: in fact, in this 
case, the system behaves like a single giant oscillator and 
becomes thus integrablc. 

If we consider again eq.(3) describing a system of cou- 
pled driven/damped pendula, one can immediately ver- 
ify that eq.(5) represents its overdampcd limit, i.e. the 
case B >> 1 [6]. In this context, the natural frequen- 
cies LUi play the role of the torque term, while C ~ K 
and M = r. This common origin of both HMF and Ku- 
ramoto models from eq.(3) seems to indicate the existence 
of a non trivial link between the two oscillators systems, 
despite the non-Hamiltonian character of the latter. Ac- 
tually, their dynamics reveals many analogies. In refs. [14] 
the authors studied a generalized version of the Kuramoto 
model which share similarities with the behavior of the 
HMF model. On the other hand, in ref. [6] we studied 
analogies in the quasi-stationary behavior, i.e. the ap- 
pearance of metastable states near the phase transition. 
In the following we will compare the stationary behavior 
of the Kuramoto model with the equilibrium regime of the 
HMF model, with particular focus on the chaotic aspects. 
In the lower panels of Fig. 2 and Fig. 3 we show the asymp- 
totic behavior of the Kuramoto order parameter r as a 
function of the control parameter K for two different dis- 
tributions of the natural frequencies gioj), respectively a 
Gaussian and a uniform one, and for large systems of os- 
cillators (from N = 20000 to TV = 50000). As predicted 
by Kuramoto analysis, in both cases we observe a phase 
transition at a critical value Kc = ^TToT' ^^^ Gaussian 
distribution has mean and variance 1. therefore from 
the normalization condition follows g(0) = l/-\/(27r) and 
Kc ~ 1.59617.... The uniform distribution is selected in 
the range w S [~2, 2] therefore the normalization condi- 
tions gives 5(0) = 1/4 and Kc = 8/tt. In order to compare 
the two cases, we plot the order parameter as a function 
of K/ Kc and for several sizes of the system. Please no- 
tice that we show decreasing values of K/Kc in order to 
better compare Kuramoto data with those of the HMF 
model. An average over 10 realizations has been consid- 
ered for each dot. One immediately recognizes a different 
kind of transition: a continuous (2nd-order-like) one, for 
the Gaussian g{Lij) (Fig. 2) and an abrupt (Ist-ordcr-like) 
one, for the uniform g{ijj) (Fig. 3). Correspondingly, two 
different behaviors of the LLE (calculated as in the pre- 
vious section following ref. [10]) were also observed: they 
are plotted in the upper panels of Figs. 2 and 3 and clearly 
show that the LLE can be considered as a good dynamical 
indicator of the phase transitions. In fact, in both cases 
we observe a pronounced peak around the transition. But, 
while it slowly decreases for K > Kc in the Gaussian g{u!) 
case, in the uniform one it goes to zero abruptly just af- 
ter the critical point. In both cases the chaotic regime, 
characterized by a positive finite LLE, seems to be related 
with the simultaneous presence of drifting and locking os- 
cillators, i.e. with the existence of partially synchronized 
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Fig. 4: Temporal evolution of the order parameter r{t) near 
the phase transition for several runs and different distributions 
g{io): the Gaussian one, panels (a), (b) and (c) and the uniform 
one, panels (d), (e) and (f). Metastable states are visible in 
panels (a) and (f). See text. 



asymptotic stationary states, and seems not to depend on 
the size of the system (as happened below the critical en- 
ergy in the HMF model). On the other hand, as expected, 
the LLE vanishes for small or high values of the coupling, 
being in those cases the system completely homogeneous 
or fully synchronized (i.e., in both the cases, integrable). 

These results confirm previous studies [7] concerning the 
dependence of the transition order on the g{ui) distribu- 
tion, and extend the investigation of Maistrenko et al. [8] 
without any contradiction. In particular, in the latter, the 
authors show that phase chaos in Kuramoto model arises 
as soon as A^ = 4 or more oscillators interact. But, even if 
they compute the entire Lyapunov spectrum, indeed they 
take into account only relatively small system sizes (up 
to A'' < 200) and do not distinguish between different 
kinds of phase transition. Furthermore, they mainly con- 
sider the so-called symmetric Kuramoto model, where the 
natural frequencies uji are symmetrically allocated around 
the mean frequency fl: the latter is a very peculiar case, 
which gives rise to a very sharp Ist-order-like phase tran- 
sition for the order parameter, with a corresponding LLE 
that is zero for all the values of the coupling except for 
a very narrow zone around the phase transition, which 
seems to vanish increasing the size of the system. Numer- 
ical results for the symmetric g{w) distribution are shown 
in the insets of Fig. 3, where the sharp transition in the 
order parameter is clearly evident (lower inset), together 
with the correspondent size-dependent LLE behavior (up- 
per inset). This distribution is however very peculiar and 
not very realistic, although easier to deal with from an 
analytical point of view. On the other hand, compared 
with those of Fig.l, the plots of Fig. 2 seem to indicate 
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Fig. 5; (Lower panel) The asymptotic order parameter r of the 
Kuramoto model is plotted as a function of the coupling K 
for a system of A*' = 20000 oscillators and for several bounded 
Gaussian distributions (?(a^), with different standard deviations 
a and lo € [—3, 3], Increasing cr (curves from right to left), the 
(?(a;)'s change continuously from a Gaussian to a uniform dis- 
tribution and, in correspondence, the phase transition changes 
continuously from a 2nd-order-like to a Ist-order-like one, while 
the critical value of the coupling Kc shift to the left. (Upper 
panel) The Largest Lyapunov Exponent is plotted as function 
of K. Each dot represents an average over 10 runs. See text. 



Fig. 6: For the Kuramoto model we show the phase diagram 
of K versus a in the case with A^ = 20000. A well defined 
critical line, changing asymptotically from 2nd-order-like to- 
wards Ist-order-like by increasing the value of a, separates the 
fully synchronized phase from the incoherent one. The par- 
tially synchronized regime (characterized by a positive largest 
Lyapunov exponent and < r < 0.8) is also indicated (shaded 
area) between the two. In the insets, the bounded giuS) distri- 
butions (with LD £ [—3,3]) used in the simulations are plotted 
for three increasing values of a. See text. 



that the Gaussian g{w) choice for the natural frequencies 
distribution yields a phase transition and a chaotic behav- 
ior qualitatively analogous to that one found in the HMF 
model. Comparing Fig. 2 and Fig. 3, it clearly appears 
that, at variance with what happens in the uniform g{Lu) 
case, where both homogeneous and synchronized states 
simultaneously appear in correspondence of the abrupt 
phase transition, the Gaussian g{uj) distribution drives the 
Kuramoto system along an HMF-like continuous transi- 
tion without coexistence of different phases. In Fig. 4 we 
present several plots which confirm this interesting fea- 
ture. For a system oi N = 20000 oscillators, we draw the 
temporal evolution of the order parameter r(t) for several 
single runs as a function of three values K/ Kc near the 
phase transition, for both Gaussian (left column) and uni- 
form (right column) g{ijj). It clearly appears that in the 
latter case (panels (e) and (f)) stationary states with high 
and low asymptotic values of r coexist, while in the former 
case (panels (a), (b) and (c)) only partially synchronized 
stationary states are visible. Such a result reinforces the 
distinction between the Ist-order-like and 2nd-order-like 
dynamical phase transitions, occurring in the Kuramoto 
model in correspondence of different g{uj) distributions, 
which seem to play a very crucial role. As already noticed 
in ref. [6], in some cases (see for example panels (a) and 
(f)) also metastable states appear, for both g(uj) distribu- 
tions, analogously to the appearence of metastable qua- 
sistationary states (QSS) in the HMF model (when the 
system starts from out-of-equilibrium initial conditions). 



A more detailed study on these states and on their chaotic 
properties, also in relationship with the violation of Cen- 
tral Limit Theorem and with the metastable regime of the 
HMF model [15], is in preparation [16]. 

We conclude this section by showing that, although it is 
believed that for strictly unimodal gioj) distributions the 
Kuramoto phase transition should be 2nd-order-like [7], 
the shape of the transition asymptotically tends to a Ist- 
order-like one when using bounded Gaussian distributions 
with an increasing standard deviation a and u> G [—3,3]. 
We have adopted this range, instead of the u) G [—2, 2] pre- 
viously used for the uniform distribution, in order to have 
a suitable Gaussian for small value of a. In such a way 
one obtains, for a — 1, the Gaussian distribution used in 
Fig. 2, which becomes larger and larger by increasing the 
standard deviation. For cr = 5, one obtains an almost uni- 
form distribution, like that used in Fig. 3 (see for example 
the insets of Fig. 6). By calculating the order parameter 
and the LLE as function of K for several values of a, we 
obtain the curves shown in Fig. 5 for N = 20000. In this 
figure a continuous transformation from a 2nd-order HMF- 
like phase transition (lower panel, on the right) towards an 
abrupt Ist-order-like one (on the left) is clearly visible, to- 
gether with an analogous change in the peak of the largest 
Lyapunov exponent (upper panel), whose maximum value 
is also related with g{uj). In correspondence, the critical 
value Kc, initially depending on a, moves from right to 
left towards a final value which does not depend on a any 
more and is very near to the theoretical value Kc = 12 /tt 
predicted for a true uniform g{Lo) with cu G [—3, 3] (which 
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in principle would strictly request u = oo). At the same 
time, the region of partial sinchronization, which is mainly 
situated after the phase transition for small values of a (see 
Fig. 2), progressively shifts before the phase transition for 
increasing values of a . approaching the Ist-ordcr-like be- 
havior shown in Fig. 3. This scenario is summarized in the 
plot of Fig. 6, where the synchronization phase diagram of 
K versus u for the Kuramoto model is shown. Please no- 
tice that this diagram is schematic and not universal since 
it depends on the range of giyi) and is likely affected also 
by unavoidable finite size effects. We report in the three 
insets examples of g{uS) distributions for a = 1,3,5. The 
2nd-order-like critical line, drawn as a full line, separates 
the incoherent phase, with vanishing values for both r and 
the LLE, from the fully synchronized one, characterized 
by a large value of the order parameter (r > 0.8) and, 
again, by a vanishing LLE. Just around the critical line 
we found the partially synchronized regime, with positive 
LLE and values < r < 0.8. As a final remark, we notice 
that in Ref. [17] a similar phase diagram was shown for the 
HMF model. In that case the authors considered the plane 
U versus Mq, being the latter a parameter which specify 
the class of out-of-cquilibrium initial conditions leading to 
metastable quasistationary states. Such a plane was sep- 
arated into two parts by a critical line, indicating both 
2nd-order and Ist-order phase transitions from an homo- 
geneous QSS regime to a magnetized one. Despite the 
different context, we think that this analogy could be con- 
sidered a further point of contact between the HMF and 
the Kuramoto scenarios. 

Conclusions. — We presented new numerical evi- 
dence of the presence of chaotic behavior in the Kuramoto 
model for very large system sizes, discussing the analogies 
with the Hamiltonian Mean Field (HMF) model. We stud- 
ied the phase transition features and the LLE behavior for 
both models. The latter can be also regarded as the dis- 
sipative and the conservative version of a more general 
model of coupled driven/damped pendula. Our simula- 
tions confirm that two different kinds of dynamical phase 
transitions occur in the Kuramoto model, depending on 
the distribution of the natural frequencies adopted as driv- 
ing term. A uniform 3(0;) gives rise to a sharp Ist-order- 
like transition, where both homogeneous and synchronized 
stationary states coexist. Instead, a Gaussian g{uj) yields 
a continuous 2nd-order-like transition, very similar to the 
true thermodynamical phase transition observed in the 
HMF model. On the other hand, the presence in the Ku- 
ramoto model of a peak observed in the LLE correspond- 
ingly to the critical region and regardless of the kind of 
distribution g{u!) reflects the fact that in this region the 
competition between locked and drifting oscillators acti- 
vates a microscopic chaotic dynamics which is a good dy- 
namical indicator of the phase transition. Again, such a 
chaotic behavior shows many analogies with the one ob- 
served in the HMF model, which exhibits as well a peak 
just before the critical point, where there are large fluctu- 



ations in the main thermodynamical quantities character- 
izing the macroscopic phase transition. 
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